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Abstract. We present an efficient algorithm for solving local linear sys¬ 
tems with a boundary condition using the Green’s function of a connected 
induced subgraph related to the system. We introduce the method of us¬ 
ing the Dirichlet heat kernel pagerank vector to approximate local solu¬ 
tions to linear systems in the graph Laplacian satisfying given boundary 
conditions over a particular subset of vertices. With an efficient algo¬ 
rithm for approximating Dirichlet heat kernel pagerank, our local linear 
solver algorithm computes an approximate local solution with multiplica¬ 
tive and additive error e by performing 0(e -5 s 3 log(s 3 e -1 ) logn) random 
walk steps, where n is the number of vertices in the full graph and s is 
the size of the local system on the induced subgraph. 
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1 Introduction 

There are a number of linear systems which model flow over vertices of a graph 
with a given boundary condition. A classical example is the case of an electrical 
network. Flow can be captured by measuring electric current between points 
in the network, and the amount that is injected and removed from the system. 
Here, the points at which voltage potential is measured can be represented by 
vertices in a graph, and edges are associated to the ease with which current 
passes between two points. The injection and extraction points can be viewed as 
the boundary of the system, and the relationship of the flow and voltage can be 
evaluated by solving a system of linear equations over the measurement points. 

Another example is a decision-making process among a network of agents. 
Each agent decides on a value, but may be influenced by the decision of other 
agents in the network. Over time, the goal is to reach consensus among all the 
agents, in which each agrees on a common value. Agents are represented by 
vertices, and each vertex has an associated value. The amount of influence an 
agent has on a fellow agent is modeled by a weighted edge between the two 

* An extended abstract appeared in Proceedings of WAW (2013) [6]. 
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representative vertices, and the communication dynamics can be modeled by a 
linear system. In this case, some special agents which make their own decisions 
can be viewed as the boundary. 

In both these cases, the linear systems are equations formulated in the graph 
Laplacian. Spectral properties of the Laplacian are closely related to reachabil¬ 
ity and the rate of diffusion across vertices in a graph [ 3 J. Laplacian systems 
have been used to concisely characterize qualities such as edge resistance and 
the influence of communication on edges [ 23 ]. There is a substantial body of 
work on efficient and nearly-linear time solvers for Laplacian linear systems 

see also HZ]). 

The focus of this paper is a localized version of a Laplacian linear solver. In 
a large network, possibly of hundreds of millions of vertices, the algorithms we 
are dealing with and the solutions we are seeking are usually of finite support. 
Here, by finite we mean the support size depends only on the requested output 
and is independent of the full size of the network. Sometimes we allow sizes up 
to a factor of log(n), where n is the size of the network. 

The setup is a graph and a boundary condition given by a vector with speci¬ 
fied limited support over the vertices. In the local setting, rather than computing 
the full solution we compute the solution over a fraction of the graph and de 
facto ignore the vertices with solution values below the multiplicative/additive 
error bound. In essence we avoid computing the entire solution by focusing com¬ 
putation on the subset itself. In this way, computation depends on the size of 
the subset, rather than the size of the full graph. We distinguish the two cases as 
“global” and “local” linear solvers, respectively. We remark that in the case the 
solution is not “local,” for example, if all values are below the error bound, our 
alogrithm will return the zero vector - a valid approximate solution according 
to our definition of approximation. 

In this paper, we show how local Laplacian linear systems with a boundary 
condition can be solved and efficiently approximated by using Dirichlet heat 
kernel pagerank, a diffusion process over an induced subgraph. We will illus¬ 
trate the connection between the Dirichlet heat kernel pagerank vector and the 
Green’s function, or the inverse of a submatrix of the Laplacian determined 
by the subset. We also demonstrate the method of approximation using ran¬ 
dom walks. Our algorithm approximates the solution to the system restricted to 
the subset S by performing O (7~ 2 e _3 s 3 log 2 (s 3 7 _1 ) log?r) random walk steps, 
where 7 is the error bound for the solver and e is the error bound for Dirich¬ 
let heat kernel pagerank approximation, and s denotes the size of S. We as¬ 
sume that performing a random walk step and drawing from a distribution 
with finite support require constant time. With this, our algorithm runs in time 
O (7 -2 e -3 log 4 (n) log 2 (7 _1 log 3 (n))) when the support size of the solution is 
O(logn). Note that in our computation, we do not intend to compute or approx¬ 
imate the matrix form of the inverse of the Laplacian. We intend to compute an 
approximate local solution which is optimal subject to the (relaxed) definition 
of approximation. 
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1.1 A Summary of the Main Results 

We give an algorithm called Local Linear SOLVERfor approximating a local 
solution of a Laplacian linear system with a boundary condition. The algorithm 
uses the connection between the inverse of the restricted Laplacian and the 
Dirichlet heat kernel of the graph for approximating the local solution with a 
sampling of Dirichlet heat kernel pagerank vectors (heat kernel pagerank re¬ 
stricted to a subset S). It is shown in Theorem [ 4 ] that the output of Local 
Linear Solver approximates the exact local solution xs with absolute error 
0(7||6|| + [Jicsll) for boundary vector b with probability at least 1 — 7. 

We present an efficient algorithm for approximating Dirichlet heat kernel 
pagerank vectors, ApproxDirHKPR. The algorithm is an extension of the algo¬ 
rithm in [ 7 ]. The definition of e-approximate vectors is given in Section [ 5 j We 
note that this notion of approximation is weaker than the classical notions of 
total variation distance among others. Nevertheless, this “relaxed” notion of ap¬ 
proximation is used in analyzing PageRank algorithms (see [ 3 j, for example) for 
massive networks. 

The full algorithm for approximating a local linear solution, GreensSolver, 
is presented in Section [G] The algorithm is an invocation of Local Linear 
Solver with the ApproxDirHKPR called as a subroutine. The full agorithm re¬ 
quires O (7 ~ 2 e~ 3 s 3 log 2 (s 3 7 -1 ) logn) random walk steps by using the algorithm 
ApproxDirHKPR with a slight modification. Our algorithm achieves sublinear time 
after preprocessing which depends on the size of the support of the boundary 
condition. The error is similar to the error of ApproxDirHKPR. 

It is worth pointing out a number of ways our methods can be generalized. 
First, we focus on unweighted graphs, though extending our results to graphs 
with edge weights follows easily with a weighted version of the Laplacian. Second, 
we require the induced subgraph on the subset S be connected. However, if the 
induced subgraph is not connected the results can be applied to components 
separately, so our requirement on connectivity can be relaxed. Finally, we restrict 
our discussion to linear systems in the graph Laplacian. However, by using a 
linear-time transformation due to m for converting a symmetric, diagonally 
dominant linear system to a Laplacian linear system, our results apply to a 
larger class of linear systems. 

1.2 Organization 

In Section [ 2 j we give definitions and basic facts for graph Laplacian and heat 
kernel. In Section [ 3 ] the problem is introduced in detail and provides the set¬ 
ting for the local solver. The algorithm, Local Linear Solver, is presented in 
Section [ 4 j After this, we extend the solver to the full approximation algorithm 
using approximate Dirichlet heat kernel pagerank. In Section [ 5 j we give the def¬ 
inition of local approximation and analyze the Dirichlet heat kernel pagerank 
approximation algorithm. In Section [6j the full algorithm for computing an ap¬ 
proximate local solution to a Laplacian linear system with a boundary condition, 
GreensSolver, is given. Finally in Section [7] we illustrate the correctness of the 
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algorithm with an example network and specified boundary condition. The ex¬ 
ample demonstrates visually what a local solution is and how GreensSolver suc¬ 
cessfully approximates the solution within the prescribed error bounds when the 
solution is sufficiently local. 

2 Basic Definitions and Facts 

Let G be a simple graph given by vertex set V = V(G) and edge set E = E(G). 
Let u ~ v denote {it, v} £ E. When considering a real vector / defined over the 
vertices of G, we say / £ R v and the support of / is denoted by supp(/) = {u £ 
V : f(v) ^ 0}. For a subset of vertices S C V, we say s = |S| is the size of S and 
use / £ to denote vectors defined over S. When considering a real matrix 
M defined over V, we say M £ M v xV , and we use Ms to denote the submatrix 
of M with rows and columns indexed by vertices in S. Namely, Ms £ R SxS . 
Similarly, for a vector / € we use fs to mean the subvector of / with entries 
indexed by vertices in S. The vertex boundary of S is S(S) = {u £ V\S : {u, u} £ 
E for some v £ S}, and the edge boundary is d(S) = {{ u, v} £ E : u £ S,v ^ S}. 

2.1 Graph Laplacians and heat kernel 

For a graph G, let A be the indicator adjacency matrix A £ {0, l} VxV for which 
A U v = 1 if and only if {u, u} £ E. The degree of a vertex v is the number of 
vertices adjacent to it, d v = |{u £ V\A uv = 1} |. Let D be the diagonal degree 
matrix with entries D vv = d v on the diagonal and zero entries elsewhere. The 
Laplacian of a graph is defined to be L = D — A. The normalized Laplacian, 
C = D-^LD- 1 / 2 , is a degree-nomalized formulation of L, given by 

(1 if u = v, 

[ 0 otherwise. 

Let P = D^^^A be the transition probability matrix for a random walk on 
the graph. Namely, if v is a neighbor of u, then P(u,v ) = 1 /d u denotes the 
probability of moving from vertex u to vertex v in a random walk step. Another 
related matrix of significance is the Laplace operator, A = I — P. We note that 
C is similar to A. 

The heat kernel of a graph is defined for real t > 0 by 

H t = e~ tc . 

Consider a similar matrix, denoted by H t = e~ tA = D~ 1 ^ 2 'H t D 1 ^ 2 . For a given 
t £ R + and a preference vector / £ , the heat kernel pagerank is defined by 

Ptj = f T H u 

where f T denotes the transpose of /. When / is a probability distribution on V, 
we can also express the heat kernel pagerank as an exponential sum of random 
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walks. Here we follow the notation for random walks so that a random walk step 
is by a right multiplication by P: 



2.2 Laplacian Linear System 

The examples of computing current flow in an electrical network and consensus 
in a network of agents typically require solving linear systems with a boundary 
condition formulated in the Laplacian L = D — A, where D is the diagonal 
matrix of vertex degrees and A is the adjacency matrix of the network. The 
problem in the global setting is the solution to Lx = b, while the solution x is 
required to satisfy the boundary condition b in the sense that x(u) = b(u) for 
every vertex v in the support of b. Because our analysis uses random walks, we 
use the normalized Laplacian C = D~ l ^ 2 LD~ 1 ^ 2 . We note that the solution x 
for Laplacian linear equations of the form Cx = b is equivalent to solving Lx = b 
if we take x = D~ 1 ^ 2 x and b = D l l 2 b. Specifically, our local solver computes 
the solution x restricted to S', denoted xs , and we do this by way of the discrete 
Green’s function. 

Example. To illustrate the local setting, we expand upon the problem of a 
network of decision-making agents. Consider a communication network of agents 
in which a certain subset of agents / C V are followers and an adjacent subset 
l C V \ f are leaders (see Figure [l]). Imagine that the decision values of each 
agent depend on neighbors as usual, but also that the values of the leaders are 
fixed and will not change. Specifically, let d v denote the degree of agent v, or the 
number of adjacent agents in the communication network, and let x be a vector 
of decision values of the agents. Suppose every follower Vf continuously adjusts 
their decision according to the protocol: 


1 ^ x ( u ) 

X(Vf) = x{Vf) -N - 



while every leader i>i remains fixed at b(vi). Then the vector of decision values x 
is the solution to the system Cx = b , where x is required to satisfy the boundary 
condition. 

In our example, we are interested in computing the decision values of the 
followers of the network where the values of the leaders are a fixed boundary 
condition, but continue to influence the decisions of the subnetwork of followers. 

3 Solving Local Laplacian Linear Systems with a 
Boundary Condition 

For a general connected, simple graph G and a subset of vertices S, consider the 
linear system Cx = 6, where the vector b has non-empty support on the vertex 
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Fig. 1: A communication network of agents where the leaders (in purple) have 
fixed decisions and the followers (in red) compute their decisions based on the 
leaders and the subnetwork of followers. The local solution would be the decisions 
of the followers. 


boundary of S. The global problem is finding a solution x that agrees with b , in 
the sense that x(v) = b(v) for every vertex b in the support of b. In this case we 
say that x satisfies the boundary condition b. 

Specifically, for a vector b £ lA, let S denote a subset of vertices in the 
complement of supp(fr). Then b can be viewed as a function defined on the 
vertex boundary 5(S) of S and we say b is a boundary condition of S. Here we 
will consider the case that the induced subgraph on S is connected. 

Definition 1. Let G be a graph and let b be a vector b £ R' / over the vertices 
of G with non-empty support. Then we say a subset of vertices S C V is a 
6 -boundable subset if 

(i) SCV\suppib), 

(ii) S(S) D suppib) ^ 0, 

(Hi) the induced subgraph on S is connected and 6{S) ^ 0. 

We note that condition ( [Til] ) is required in our analysis later, although the 
general problem of finding a local solution over S can be dealt with by solving 
the problem on each connected component of the induced subgraph on S in¬ 
dividually. We remark that in this setup, we do not place any condition on b 
beyond having non-empty support. The entries in b may be positive or negative. 

The global solution to the system Cx = b satisfying the boundary condition 
b is a vector x £ with 


x{v) 



if v £ S 
iiv^S 


( 1 ) 
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for a 6 -boundable subset S. The problem of interest is computing the local so¬ 
lution for the restriction of x to the subset S, denoted xs- 

The eigenvalues of Cs are called Dirichlet eigenvalues, denoted Ai < A 2 < 
• • • < A s where s = |5|. It is easy to check (see jlj) that 0 < Ai < 2 since we 
assume S(S) 0. Thus Cg 1 exists and is well defined. In fact, s~ 3 < Ai < 1. 

Let As,ss be the s x |<5(5)| matrix by restricting the columns of A to S(S) and 
rows to S. Requiring S' to be a 6 -boundable subset ensures that the inverse C g 1 
exists [4j. Then the local solution is described exactly in the following theorem. 

Theorem 1. In a graph G, suppose b is a nontrivial vector in KA and S is a 
b-boundable subset. Then the local solution to the linear system Cx = b satisfying 
the boundary condition b satisfies 

x s = ^s 1 (-C ) s 1 ^ 2 As,ssD S g/~bss). (2) 


Proof. The vector b\ := Dg 1 ^ 2 A$sD S g^ 2 bss is defined over the vertices of S, 
and giveover the vertices of S by 


h( v ) = E 


u£ 5 (S) ,u~v 

Also, the vector jCsXs is given by, for v £ S, 


b(u) 

\/ d v d u 


C s x s (v) = x(v) - ^2 


u£S,u~v 


x(u) 
yj dydy 


By 0 and 0 , we have 


x s(v) = 


x(u) 


E 


b(u) 


c-a „ y/d v d u ., c , y/dydn 


and combining (l3|) and (|4j) , we have that xs = Cg 2 bi. 


(3) 


(4) 


□ 


3.1 Solving the local system with Green’s function 

For the remainder of this paper we are concerned with the local solution xs- 
We focus our discussion on the restricted space using the assumptions that the 
induced subgraph on S is connected and that 6(S) 0. In particular, we consider 

the Dirichlet heat kernel , which is the heat kernel pagerank restricted to S. 

The Dirichlet heat kernel is written by Hs,t and is defined as Hs,t = e~ tCs . 
It is the symmetric version of Hs,t, where Hs t t = e~ tAs = Dg 1 ^~'Hs,tD 1 J 2 . 

The spectral decomposition of Cs is 

S 

£s = ^2 

2=1 
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where P-, are the projections to the ?’th orthonormal eigenvectors. The Dirichlet 
heat kernel can be expressed as 


n s ,t = ' 


-tAo 


Let Q denote the inverse of £§■ Namely, QCs = CsQ = Is- Then 


1 




2 = 1 


From ([5]) , we see that 


2^"^’ 


(5) 


( 6 ) 


where || • || denotes the spectral norm. We call Q the Greens function , and Q can 
be related to TLs,t as follows: 

Lemma 1. Let Q be the Green’s function of a connected induced subgraph on 
S C V with s = |Sj. Let Lis t be the Dirichlet heat kernel with respect to S. Then 


pOO 

G = ^s,t dt- 

Jo 


/ 0 

Proof. By our definition of the heat kernel, 

nOO . s 


r°° r°°/ \ 

/ Hs,tdt= (£y tA *Pi) 

Jo Jo v i=1 7 

-tc 


d t 


, —tXi 


dt Pj 


1 

2=1 

= g. 


□ 


Equipped with the Green’s function, the solution ([ 2 ]) can be expressed in 
terms of the Dirichlet heat kernel. As a corollary to Theorem [l] we have the 
following. 

Corollary 1 . In a graph G, suppose b is a nontrivial vector in JR' 7 and S is a 
b-boundable subset. Then the local solution to the linear system. Cx = b satisfying 
the boundary condition b can be written as 

pOO 

xs = / Us,tbi dt, (7) 

Jo 

where b x = Dg 1/2 AsjsDffpb&s■ 

The computation of b\ takes time proportional to the size of the edge bound¬ 
ary. 
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4 A Local Linear Solver Algorithm with Heat Kernel 
Pagerank 

In the previous section, we saw how the local solution x$ to the system satisfying 
the boundary condition b can be expressed in terms of integrals of Dirichlet heat 
kernel in ([7]). In this section, we will show how these integrals can be well- 
approximated by sampling a finite number of values of Dirichlet heat kernel 
(Theorem [2]) and Dirichlet heat kernel pagerank (Corollary [ 2 ]). All norms || • || 
in this section are the L 2 norm. 

Theorem 2. Let G be a graph and C denote the normalized Laplacian of G. 
Let b be a nontrivial vector b £ lA and S a b-boundable subset, and let b\ = 
Dg 1 ^ 2 As t ssD s ^ 2 bss■ Then the local solution xs to the linear system Cx = b 
satisfying the boundary condition b can be computed by sampling TLs,tbi for r = 
7 ~ 2 log(s 7 _1 ) values. If xs is the output of this process, the result has error 
bounded by 

Iks - £s|| = 0 ( 7(1 |6iII + I|a:s|D) 

with probability at least 1 — 7 . 

We prove Theorem [ 2 ] in two steps. First, we show how the integral Q can 
be expressed as a finite Riemann sum without incurring much loss of accuracy 
in Lemma [2] Second, we show in Lemma [3] how this finite sum can be well- 
approximated by its expected value using a concentration inequality. 

Lemma 2. Let xs be the local solution to the linear system Cx = b satisfying the 
boundary condition b given in 0 - Then , for T = s 3 log(s 3 7 3 ) and N = T/'y, 
the error incurred by taking a right Riemann sum is 

N T 

IAS -J2 H S,jT/NJ^b i|| < 7(||&i|| + |ks||), 

3=1 

where 61 = Df, 1/2 A s ,ssDJg /2 bss■ 

Proof. First, we see that: 


iiHs,di = nE e_tAip ‘ii 

i 

<e-^|l£ P ‘H 

= e" tAl ( 8 ) 


where A,; are Dirichlet eigenvalues for the induced subgraph S. So the error 
incurred by taking a definite integral up to t = T to approximate the inverse is 
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the difference 


nT /»oo 

Iks - / Hs,tbi di|| = || / Hs,tbi dt|| 

Jo Jt 

/ oo 

e _tAl ||&i|| d t 

<fe- TAl ||k||. 

Then by the assumption on T the error is bounded by ||xs — Hsjbi d£|| < 

Tll&ill* 

Next, we approximate the definite integral in [0, T] by discretizing it. That is, 
for a given 7 , we choose N = T/ 7 and divide the interval [0,T] into N intervals 
of size T/N. Then a finite Riemann sum is close to the definite integral: 

II / W-sjbi dt~y^ W-sjT/Nbi 77 11 < 7 I| [ Us,tbi dt|| 

J 0 j=1 JV J 0 

< 7lks||- 

This gives a total error bounded by 7 (||£>i11 + |ks||)- □ 

N 

Lemma 3. The sum H S,jT/Nbijj can be approximated by sampling r = 

i=l 

7 ~ 2 log(s 7 _1 ) values ofTlsjT/Nbi where j is drawn from [1, TV]. With probability 
at least 1 — 7 , the result has multiplicative error at most 7 . 

A main tool in our proof of Lemma [3] is the following matrix concentration 
inequality (see [5], also variations in [1], 1, ED], EE], ESI). 


Theorem 3. Let X\, X 2 , ■ ■ ■, X m be independent random n x n Hermitian ma¬ 
trices. Moreover, assume that ||X; — E(Xj)|| < M for all i, and put v 2 = 
|| var {Xi) ||. Let X = Xi. Then for any a > 0, 

Pr(||X -E(X)|| > a) < 2„exp (■- — 2 + ° M<i/3 ) , 

where || • || denotes the spectral norm. 

Proof of Lemma [3] Suppose without loss of generality that || 6 i|| = 1. Let Y be 
a random variable that takes on the vector TLs,jT/Nbi for every j £ [1, TV] with 

N r 

probability 1/N. Then E(Y) = )C Tls,jT/Nb-\ ■ Let X = J2 Xj where each Xj 

j= 1 ’ i —1 

is a copy of Y, so that E(X) = rE(Y). 

Now consider Y to be the random variable that takes on the projection matrix 
'H-S,jT/Nbi('Hs,jT/Nbi) T for every j £ [1, AT] with probability 1/N, and X is the 
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sum of r copies of Y. Then we evaluate the expected value and variance of X as 
follows: 


||E(X)||=r||E(Y)|| 

N 

11Var(X)11 =r||Var(Y)|| < \\—'^2'Hs,jT/Nbi('Hs,jT/Nbi) T \\'Hs,jT/Nbi\\ 2 

3=1 

<r ||E(Y)||. 


We now apply Theorem [3] to X. We have 


Pr(||X - E(X)11 > 7||E(X)||) < 2sexp 


< 2s exp 


< 2s exp 


7 2 ||E(X)|| 2 \ 

2Var(X) + ) 

7 2 r 2 ||E(Y)|| \ 
r + 2'yrM/3 ) 


Therefore we have Pr(||X—E(X)| | > y| |E(X)| |) < 7 if we choose r > 7 2 log(s7 - 1 ). 
Further, this implies the looser bound: 

Pr(||X-E(X)|| >7||E(X)||) <7. 

Then E(Y) = yE(X) is close to y X and 


N 1 1 N 1 
I n SJT/Nbl — - -X\ | < 7 | | HsjT/Nb 1 

3=1 3=1 


N T T N 

\^2'H S , J T/Nbl^ r - —^11 < l\\^2T-is,jT/Nbi 


3=1 


3=1 


N' 

T 
N 1 


with probability at least 1 — 7, as claimed. □ 

Proof of Theorem Q Let X be the sum of r samples of TLs jT/Nbi with j drawn 
from [0, AT], and let xs = yX. Then combining Lemmas [2] and [3J we have 

N T 

Iks - X S \\ < 7(1 kill + Iks]] + \\J2 n S,3T/Nbl — \\) 

3=1 

<o( 7 (INI + lksll)). 

By Lemma [3j this bound holds with probability at least 1 — 7 . □ 

The above analysis allows us to approximate the solution xs by sampling 
T~Ls,tb\ for various t. The following corollary is similar to Theorem [2] except we 
use the asymmetric version of the Dirichlet heat kernel which we will need later 
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for using random walks. In particular, we use Dirichlet heat kernel pagerank 
vectors. Dirichlet heat kernel pagerank is also defined in terms of a subset S 
whose induced subgraph is connected, and a vector / £ R 5 by the following: 

PS,tj = ( 9 ) 


Corollary 2. Let G be a graph and C denote the normalized Laplacian of G. 
Let b be a nontrivial vector b £ R y and S be a b-boundable subset. Let b 2 = 
(Dg 1 ^ 2 As, 5 sD S g^ 2 bss) T D 1 ^ 2 . Then the local solution xs to the linear system 
Cx = b satisfying the boundary condition b can be computed by sampling ps,t,b 2 
for r = 7 ^ 2 log(s 7 ^ 1 ) values. If xs is the output of this process, the result has 
error bounded by 

Iks -zsll = 0(7(1 l^i II + Iksll)), 

where b\ = D s ' AssD ss bss, with probability at least 1 — 7 . 

Proof. First, we show how x$ can be given in terms of Dirichlet heat kernel 
pagerank. 


roc 

Xs = / b{Us,t 
Jo 

= r^Hs,D^) d t 

Jo 

rOO 

= / b 2 Hs,tL>g 1 ^ 2 dt, where b 2 = b{D l J 2 
Jo 


PS,t,bo 


df D 


- 1/2 


and we have an expression similar to ([t]) . Then by Lemma [ 2 J x"g is close to 
N 

S Ps,jT/N,b 2 1 ^ 2 with error bounded by 0(7(1 |&i|| + | ks||)) ■ From Lemma 3 
3 =1 

this can be approximated to within 0 ( 7 ||a:s||) multiplicative error using r = 
7 -2 log(s 7 _1 ) samples with probability at least 1 — 7 . This gives total additive 
and multiplicative error within 0( 7 ). □ 


4.1 The Local Linear Solver Algorithm 

We present an algorithm for computing a local solution to a Laplacian linear 
system with a boundary condition. 

Theorem 4. Let G be a graph and C denote the normalized Laplacian of G. 
Let b be a nontrivial vector b £ R v , S a b-boundable subset, and let b\ = 
Dg 1 ^ 2 As,5sDss^ 2 bss■ For the linear system Cx = b, the solution x is required 
to satisfy the boundary condition b, and let xs be the local solution. Then the 
approximate solution x output by the Local Linear Solver algorithm has an 
error bounded by 

Iks — x 11 = o( 7 (INI + Iksll)) 

with probability at least 1 — 7 . 
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Algorithm 1 Local Linear Solver 

input: graph G, boundary vector b € R v , subset S C V, solver error parameter 
0 < 7 < 1 . 

output: an approximate local solution x with additive and multiplicative error 7 to 
the local system xs = Qbi satisyfing the boundary condition b. 

1: a<-\S\ 

2 : initialize a 0 -vector x of dimension s 
3: 61 «- D- 1/2 A s , SS D^ /2 bss 
4: fe 2 <- bf D 1 / 2 
5: T <— s 3 log(s 3 7 _1 ) 

6 : N <— T /7 

7: r 7 -2 log(s 7 _1 ) 

8 : for i = 1 to r do 

9: draw j from [1, TV] uniformly at random 

10: Xi <— pS,jT/N,b 2 

11 : x «— x + Xi 

12: end for 

13: return T/r ■ xDg 1 ^ 2 


Proof. The correctness of the algorithm follows from Corollary [2] □ 

The algorithm involves r = 7 _2 log(s 7 ” 1 ) Dirichlet heat kernel pagerank 
computations, so the running time is proportional to the time for computing 
b 2 e~ TAs for T = s 3 log(s 3 7 _1 ). 

In the next sections, we discuss an efficient way to approximate a Dirichlet 
heat kernel pagerank vector and the resulting algorithm GreensSolver that 
returns approximate local solutions in sublinear time. 

5 Dirichlet Heat Kernel Pagerank Approximation 
Algorithm 

The definition of Dirichlet heat kernel pagerank in is given in terms of a 
subset S and a vector / £ Our goal is to express this vector as the station¬ 
ary distribution of random walks on the graph in order to design an efficient 
approximation algorithm. 

Dirichlet heat kernel pagerank is defined over the vertices of a subset S as 
follows: 


Ps,t,f = f T H s , t = f T e~ tAs 

OO . h 

= T, e ~ t yf Tp s- 

k—0 


= / J e 


T p -t(I s -Ps) 


That is, it is defined in terms of the transition probability matrix P$ - the restric¬ 
tion of P where P describes a random walk on the graph. We can interpret the 
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matrix Pg as the transition probability matrix of the following so-called Dirichlet 
random walk: Move from a vertex u in S to a neighbor v with probability l/d u . 
If v is not in S, abort the walk and ignore any probability movement. Since we 
only consider the diffusion of probability within the subset, any random walks 
which leave S cannot be allowed to return any probability to S. To prevent this, 
random walks that do not remain in S are ignored. 

We recall some facts about random walks. First, if g is a probabilistic function 
over the vertices of G, then g T P k is the probability distribution over the vertices 
after performing k random walk steps according to P starting from vertices 
drawn from g. Similarly, when / is a probabilistic fuction over S , f T Pg is the 
distribution after k Dirichlet random walk steps. Consider a Dirichlet random 
walk process in which the number of steps taken, k (where steps are taken 
according to a Dirichlet random walk as described above), is a Poisson random 
variable with mean t. That is, k steps are taken with probability pt.{k) = e~ l -^. 
Then, the Dirichlet heat kernel pagerank is the expected distribution of this 
process. 

In order to use random walks for approximating Dirichlet heat kernel pager¬ 
ank, we perform some preprocessing for general vectors f £ R s . Namely, we 
do separate computations for the positive and negative parts of the vector, and 
normalize each part to be a probability distribution. 

Given a graph and a vector / £ M 5 / the algorithm ApproxDirHKPR computes 
vectors that e-approximate the Dirichlet heat kernel pagerank ps.t.j satisfying 
the following criteria: 

Definition 2. Let G be a graph and let S C V be a subset of vertices. Let 
f £ be a probability distribution vector over the vertices of S and let ps,tj be 
the Dirichlet heat kernel pagerank vector according to S, t, and f. Then we say 
that v £ M. s is an e-approximate Dirichlet heat kernel pagerank vector if 

1. for every vertex v in the support of v, \ps,t,f(v) — v(v)\ < e • ps,tj( v ), an d 

2. for every vertex with v(v) = 0, it must be that ps,tj( v ) < e - 

When f is a general vector, an e-approximate Dirichlet heat kernel pagerank 
vector has an additional additive error o/e||/||i by scaling, where || • ||i denotes 
the Li norm. 

For example, the zero-vector is an e-approximate of any vector with all entries 
of value < e. We remark that for a vector / with L\ norm 1, the Dirichlet heat 
kernel pagerank vector ps,tj has at most 1/e entries with values at least e. Thus 
a vector that e-approximates ps,tj has support of size at most 1 /e. 

The time complexity of ApproxDirHKPR is given in terms of random walk 
steps. As such, the analysis assumes access to constant-time queries returning 
(i) the destination of a random walk step, and (ii) a sample from a distribution. 

Theorem 5. Let G be a graph and S a proper vertex subset such that the induced 
subgraph on S is connected. Let f be a vector f £ t £ R + , and 0 < e < 
1. Then the algorithm ApproxDirHKPR (G,t, /, S,e) outputs an e-approximate 
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Algorithm 2 ApproxDirHKPR (G,t, f, S,e) 

input: a graph G, t £ R + , vector / £ R s , subset S C V, error parameter 0 < e < 1. 
output: p, an e-approximation of ps,t,f- 

1: s i— |S| 

2: initialize 0-vector p of dimension s 
3: /+ «— the positive portion of / 

4: /_ <— the negative portion of / so that / = /+ — /_ 

5: f' + t— /+/||/+||i > normalize /+ to fee a probability distribution vector 

6: /i <— /_/||/_||i > normalize /_ to 6e a probability distribution vector 

7: r <- if logn 
8: for r iterations do 

9: choose a starting vertex u\ according to the distribution vector 

10: fc ~ Poiss(t) > choose A: with probability e -4 ^ 

11: fc <— min{k, t/e} 

12: simulate k steps of a P = D^ 1 A random walk 

13: if the random walk leaves S then: 

14: do nothing for the rest of this iteration 

15: else 

16: let Vi be the last vertex visited in the walk 

17: p[vi] <- p[vi] + 11 /+11 1 

18: end if 

19: choose a starting vertex U 2 according to the distribution vector f'_ 

20: k ~ Poiss(t) > choose k with probability 

21: k min{k, t/e} 

22: simulate k steps ofaP = D~ 1 H random walk 

23: if the random walk leaves S then: 

24: do nothing for the rest of this iteration 

25: else 

26: let V 2 be the last vertex visited in the walk 

27: p[v 2 \ <r- p[v 2 ] + ||fl—||i 

28: end if 

29: end for 
30: 

31: p <— p/r 

32: return p 
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Dirichlet heat kernel pagerank vector ps,tj with probability at least 1 — e. The 
running time of ApproxDirHKPR is 0^e _4 flogn^, where the constant hidden in 
the big-0 notation reflects the time to perform a random walk step. 

Our analysis relies on the usual ChernofT bounds restated below. They will 
be applied in a similar fashion as in [3]. 

Lemma 4 ([3]). Let Xi be independent Bernoulli random variables with X = 
Xi. Then, 

i -1 

1. for 0 < e < 1, Pr(A' < (1 — e)rE(A)) < exp(—/^rE(A)) 

2. for 0 < e < 1, Pr(A > (1 + e)rE(X)) < exp(— ^rE(A)) 

3. for c> 1, Pr(A' > (1 + c)?’E(A)) < exp(— |rE(X)). 

Proof of Theorem [5| For the sake of simplicity, we provide analysis for the pos¬ 
itive part of the vector, / := / + , noting that it is easily applied similarly to the 
negative part as well. 

The vector f = //||/||i is a probability distribution and the heat kernel 
pagerank p'g t f = As,t,//||/||i can be interpreted as a series of Dirichlet random 

walks in which, with probability f ,T Pg is contributed to p'stf - This is 

demonstrated by examining the coefficients of the terms, since 


The probability of taking k ~ Pois[t ) steps such that k >t/e is less than e by 
Markov’s inequality. Therefore, enforcing an upper bound of K = t/e for the 
number of random walk steps taken is enough mixing time with probability at 
least 1 — e. 

For k < t/e, our algorithm approximates f' T Pg by simulating k random walk 
steps according to P as long as the random walk remains in S. If the random 
walk ever leaves S, it is ignored. To be specific, let Xf be the indicator random 
variable defined by Xf = 1 if a random walk beginning from a vertex u drawn 
from /' = //||/|1 1 ends at vertex v in k steps without leaving S. Let A'” be 
the random variable that considers the random walk process ending at vertex v 
in at most k steps without leaving S. That is, X v assumes the vector Xf with 
probability Namely, we consider the combined random walk 

A " = E e" ’PXI 

k<t/e 

Now, let p{k)s,t,f be the contribution to the heat kernel pagerank vector 
p' s t f of walks of length at most k. The expectation of each X v is p{k)s,t,f{v)- 
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Then, by Lemma [4j 

Pi(X v < (1 - e)p(k) s ,tj(v) ■ r) < exp(-p(k) s ,tj(v)re 2 /2) 

= exp(-(8 /e)p(k) s ,tj(v) log n) 

<n- 4 

for every component with p' St j(v) > e, since then p(k)s t t,f(v) > e/2. Similarly, 

Pi(X v > (1 + e)p(k) s ,tj(v) ■ r) < exp(-p(k) s ,t,f{v)re 2 /4:) 

= exp(-(4/e)p(fc) 5itj/ (u) log n) 

< n~ 2 . 

We conclude the analysis for the support of p' St ^ by noting that Ps,t.,f = \ X v , 
and we achieve an e-multiplicative error bound for every vertex v with p' s t j(v) > 
e with probability at least 1 — 0(n~ 2 ). 

On the other hand, if p' s t f{v) < e, by the third part of Lemma|4j Pr(ps jt j(v) > 
2e) < n~ s ^ 2 . We conclude that, with high probability, ps,tj{v) < 2e. 

Finally, when / is not a probability distribution, the above applies to f = 
//I I/ 111 - Let p's t f be the output of the algorithm using /' = //||/||i and p' s t j 
be the corresponding Dirichlet heat kernel pagerank vector Ps,t,p - The full error 
of the Dirichlet heat kernel pagerank returned is 

\\Ps,tj ~ Ps,tj\\i < Illl/lliPs.t,/ - ll/lliPs,t,/lli 
< \\f\\l\\p's,tj ~ p's,t,/111 

<e||/lli||pWHi 

= e||/||i- 

For the running time, we use the assumptions that performing a random 
walk step and drawing from a distribution with finite support require constant 
time. These are incorporated in the random walk simulation, which dominates 
the computation. Therefore, for each of the r rounds, at most K steps of the 
random walk are simulated, giving a total of rK = O^iflogn • t/ej = 0(t ) 
queries. □ 

6 The GreensSolver Algorithm 

Here we present the main algorithm, GreensSolver, for computing a solution to 
a Laplacian linear system with a boundary condition. It is the Local Linear 
Solver algorithmic framework combined with the scheme for approximating 
Dirichlet heat kernel pagerank. The scheme is an optimized version of the algo¬ 
rithm ApproxDirHKPR with a slight modification. We call the optimized version 
SolverApproxDirHKPR. 
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Definition 3. Define SolverApproxDirHKPR (G,t, /, S, e) to be the algorithm 


after drawing k ~ Poiss(t): 
k min{k, 2 1}. 

Namely, this modification limits the length of random walk steps to at most 2 1. 

Theorem 6. Let G be a graph and S a subset of size s. Let T = s 3 log(s 3 7 _1 ), 
and let N = T/ 7 for some 0 < 7 < 1. Suppose j is a random variable drawn 
from [1, [N\] uniformly at random and let t = jT/N. Then if e > 7 , the al¬ 
gorithm SolverApproxDirHKPR returns a vector that e-approximates ps,t,f with 
probability at least 1 — e. Using the same query assumptions as Theorem [dj the 
running time of SolverApproxDirHKPR is O (w 3 flogn). 

We will use the following Chernoff bound for Poisson random variables. 

Lemma 5 ( [19] ). Let X be a Poisson random variable with parameter t. Then, 
if x > t, 

Pr(X > x) < e z-*-* lo g(*/t). 

Proof of Theorem [3J Let k be a Poisson random variable with parameter t. Sim¬ 
ilar to the proof of Theorem [5] we use Lemma [5] to reason that 

Pr(fc > 2 1) < e 2 t -*- 2 “°s( 2t /‘) 

_ g*( 1-2 log 2) 

<e, 


ApproxDirHKPR(G,t, f,S,e) with the following modification to lines \11\ and 21 


as long as t > 


Let E be the event that t < 7357 ^- The probability of E is 


Pr jT/N < 


l°g(e ) 

1 - 21 og 2 


= Pr j < 


log(e ) 

7 ( 1-2 log 2 ) 

l°g( e 1 ) 


(1 — 2 log 2 )s 3 log(s 3 7 “ 


which is less than e as long as e > (t -)*' 1 21 og “P‘ s . This holds when e > 7 . 

As before, the algorithm consists of r rounds of random walk simulation, 
where each walk is at most 2 1. The algorithm therefore makes r-2t = e -3 32f logn 
queries, requiring O (e -3 f logn) time. □ 


Below we give the algorithm GreensSolver. The algorithm is identical to 
Local Linear Solver with the exception of line [ 10 J where we use the ap¬ 
proximation algorithm SolverApproxDirHKPR for Dirichlet heat kernel pagerank 
computation. 
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Algorithm 3 GreensSolver (G, b, S, 7 , e) 

input: graph G, boundary vector b £ R v , subset S C V, solver error parameter 
0 < 7 < 1, Dirichlet heat kernel pagerank error parameter 0 < e < 1. 
output: an approximate local solution x to the local system xs = Qbi satisyfing the 
boundary condition b. 

1: s «— |S| 

2 : initialize a 0 -vector x of dimension s 
3: 61 «- D- 1/2 A s , SS D^ /2 bss 
4: fe 2 <- bf D 1 / 2 
5: T <— s 3 log(s 3 7 _1 ) 

6 : N <— T /7 

7: r ■£- 7 _ 2 log(s 7 _1 ) 

8 : for i = 1 to r do 

9: draw j from [1, TV] uniformly at random 

10: Xi <— SolverApproxDirHKPR(G, jT/N, & 2 , S, t) 

11 : x «— x + 07 

12: end for 

13: return T/r ■ xDg 1 ^ 2 


Theorem 7. Let G be a graph and C denote the normalized Laplacian of G. 
Let b be a nontrivial vector b £ R^ and S a b-boundable subset, and let b± = 
D s ' As^sDss bgs- For the linear system Lx = b, the solution x is required 
to satisfy the boundary condition b, and let xg be the local solution. Then the 
approximate solution x output by the algorithm GreensSolver satisfies the fol¬ 
lowing: 

(i) The error of x is ||xs — x|| = 0(7(||i>i|| + ||xs||) + e||& 2 11 1 ) with probability 
at least 1^7, 

(ii) The running time of GreensSolver is O ( 7 -2 e -3 s 3 log 2 (s 3 7 _1 ) logn) where 
the big-0 constant reflects the time to perform a random walk step, plus ad¬ 
ditional preprocessing time 0(|i9(S')|), where d{S) denotes the edge bound¬ 
ary of S. 

Proof. The error of the algorithm using true Dirichlet heat kernel pagerank 
vectors is 0(7(||&i|| + ||xs||)) by Corollary^ so to prove ^ we address the ad¬ 
ditional error of vectors output by the approximation of SolverApproxDirHKPR. 
By Theorem [6j SolverApproxDirHKPR outputs an e-approximate Dirichlet heat 
kernel pagerank vector with probability at least 1 — e. Let ps,t,f be the output of 
an arbitrary run of SolverApproxDirHKPR(G, t, /, S, e) . Then \\ps,tj ~Ps,t,f\\ < 
e(||ps,t,/'||i + ||/||i) = e11./11 1 by the definition of e-approximate Dirichlet heat 
kernel pagerank vectors, where f = //||/||i is the normalized vector /. This 
means that the total error of GreensSolver is 

ll^s-xH < O (7(| |i»i 11 + ||x s |D) +e||6 2 ||i. 

Next we prove |h|. The algorithm makes r = 7 -2 log(s 7 _1 ) sequential calls to 
SolverApproxDirHKPR. The maximum possible value of t is T = s 3 log(s 3 7 _1 ), 
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so any call to SolverApproxDirHKPR is bounded by O (e 3 s 3 log(s 3 j 1 )log n). 
Thus, the total running time is O ( 7 - 2 e - 3 s 3 log 2 (s 3 7 _1 ) log re). 

The additional preprocessing time of 0(|9(S)|) is for computing the vectors 
61 and 62 ; these may be computed as a preliminary procedure. □ 

We note that the running time above is a sequential running time attained 
by calling SolverApproxDirHKPR r times. However, by calling these in r parallel 
processes, the algorithm has a parallel running time which is simply the same 
as that for SolverApproxDirHKPR. 


6.1 Restricted Range for Approximation 


Since SolverApproxDirHKPR only promises approximate values for vertices whose 
true Dirichlet heat kernel pagerank vector values are greater than e, the 
GreensSolver algorithm can be optimized even further by preempting when this 
is the case. 

Figure [2] illustrates how vector values drop as t gets large. The network is the 
same example network given in Section |2.2| and is further examined in the next 
section. We let t range from 1 to T = s 3 log(s 3 7 _1 ) ~ 108739 for 7 = 0.01 and 
compute Dirichlet heat kernel pagerank vectors ps,t,f■ The figure plots L 1 norms 
of the vectors as a solid line, and the absolute value of the maximum entry in 
the vector as a dashed line. In this example, no vector entry is larger than 0.01 
for t as small as 250. 

Suppose it is possible to know ahead of time whether a vector ps,t.j will have 
negligably small values for some value t. Then we could skip the computation of 
this vector and simply treat it as a vector of all zeros. 

From ([ 8 ]), the norm of Dirichlet heat kernel pagerank vectors are mono¬ 
tone decreasing. Then it is enough to choose a threshold value if beyond which 
I \ps,t',f 111 < e, since any e-approximation will return all zeros, and treat this 
as a cutoff for actually executing the algorithm. An optimization heuristic is to 
only compute SolverApproxDirHKPR (G, t, /, S, e) if t is less than this threshold 


value if . Otherwise we can add zeros (or do nothing). That is, replace line 10 
GreensSolver with the following: 
if jT /N < t ' then 

Xi <— SolverApproxDirHKPR(G, jT/N, & 2 , S, e) 
else 

do nothing 

end if 

From @, a conservative choice for if is log(e _1 ). 


7 An Example Illustrating the Algorithm 

We return to our example to illustrate a run of the Green’s solver algorithm for 
computing local linear solutions. The network is a small communication network 

of dolphins m- 
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0 < t < T 


Fig. 2: How support values of a Dirichlet heat kernel pagerank 
vector change for different values of 1 < t < T = s 3 log(s 3 7 _1 ). 
The solid line is the L\ norm - the sum of all the support values 
- and the dashed line is the absolute value of the maximum entry 
in the vector. Note the x-axis is log-scale. 
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In this example, the subset has a good cluster, which makes it a good candi¬ 
date for an algorithm in which computations are localized. Namely, it is ideal for 
SolverApproxDirHKPR, which promises good approximation for vertices that ex¬ 
ceed a certain support threshold in terms of the error parameter e. The support 
of the vector b is limited to the set of leaders, which is the vertex boundary of 
the subset of followers, l = 6(f). The vector is plotted over the agents (vertices) 
in Figure [3] 



Fig. 3: The values of the boundary vector plotted against the 
agent IDs given in Figure [T] 


Figure [4j plots the vector values of the heat kernel pagerank vector p t ^ 
the full set of agents. Here, we use b' 2 , the n-dimensional vector: 


over 


b 2 0 ) 


62 (v) if v £ S, 
0 otherwise, 


and t = 50.0. The components with largest absolute value are concentrated in 
the subset of followers over which we compute the local solution. This indicates 
that an output of SolverApproxDirHKPR will capture these values well. 


7.1 Approximate solutions 

In the following figures, we plot the results of calls to our approximation al¬ 
gorithms against the exact solution xg using the boundary vector of Figure [3j 
The solution xg is computed by Theorem [l] and the appromimations are sample 
outputs of Local Linear SOLVERand GreensSolver, respectively. The exact 
values of xg are represented by circles, and the approximate values by triangles 
in each case. Note that we permute the indices of the vertices in the solutions 
so that vector values in the exact solution, xg are decreasing, for reading easeQ 

1 The results of these experiments as well as the source code are archived at 

http://cseweb.ucsd.edu/~osimpson/localsolverexample.html 
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Fig. 4: The node values of the full example communication net¬ 
work over a sample heat kernel pagerank vector. The red bars 
correspond to the network of followers, the purple to the leaders, 
and the white to the rest of the network. 


The result of a sample call to Local Linear Solver with error parame¬ 
ter 7 = 0.01 is plotted in Figure [5| The total relative error of this solution is 
Iplsll = 0-02, and the absolute error ||a?s — xg|| is within the error bounds 
given in Theorem[4j That is, ||xs-£s|| < 7(]|&i|| + ||xs|| + ||aVie||), where x rie 
is the solution obtained by computing the full Riemann sum (as in Lemma [2]). 

The result of a sample call to GreensSolver with parameters 7 = 0.01, e = 
0.1 is plotted in Figure[6] In this case the relative error is ~ 2.05, but the absolute 
error meets the error bounds promised in Theorem [7] point ([ 1 ]). Specifically, 

||zs - zsil < (7GI&1II + ||*s|| + \\Xrie\\) + e|| 6 2 ||i) • 


General remarks. While we have focused our analysis on solving local linear 
systems with the normalized Laplacian C as the coefficient matrix, our methods 
can be extended to solve local linear systems expressed in terms of the Laplacian 
L as well. There are numerous applications involving solving such linear systems. 
Some examples are discussed in [ 6 ], and include computing effective resistance 
in electrical networks, computing maximum flow by interior point methods, de¬ 
scribing the motion of coupled oscillators, and computing state in a network 
of communicating agents. In addition, we expect the method of approximating 
Dirichlet heat kernel pagerank in its own right to be useful in a variety of related 
applications. 


Acknowledgements. The authors would like to thank the anonymous review¬ 
ers for their comments and suggestions. Their input has been immensely helpful 
in improving the presentation of the results and clarifying details of the algo¬ 
rithm. 
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Fig. 5: The results of a run of Local Linear Solver. Two 
vectors are plotted over IDs of agents in the subset. The circles 
are exact values of a;g, while the triangles are the approximate 
values returned by Local Linear Solver. 



Fig. 6 : The results of a run of GreensSolver with 7 = 0.01, e 
0 . 1 . 
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